#include <matriz_esp.h>

#include <iomanip>
#include <algorithm>
#include <cassert>
#include <fstream>

using namespace std;

extern bool mostrar_debug;

bool primeraComponenteMenor(pair<unsigned int, double> p1, pair<unsigned int, double> p2)
{
        return p1.first < p2.first;
}

/** <CONSTRUCTORES> **/

matriz_esp::matriz_esp()
{}

/* Este contructor genera tanto la matriz M en base a la matriz pasada como parametro, a MtxM en base a M y al termino independiente B */

matriz_esp::matriz_esp(matriz<tripla_dbl> & m)
{              
        _construir_m(m);
        _construir_mtxm();
        _construir_b(m);        // Contruye el termino independiente y el vector de incognitas
       
        if(mostrar_debug)
        {
                cout << "M:    " << _m_filas << " * " << _m_colum << endl;
                cout << "MtxM: " << _mtxm_filas << " * " << _mtxm_columnas << endl;
                cout << "MtxB: " << b.size() << endl;
                cout << "X:    " << x.size() << endl << endl;
        }
}

/** <METODOS PRIVADOS> **/

void matriz_esp::_construir_m(matriz<tripla_dbl> & m)
{
        _filas_original = m.cant_filas();
        _columnas_original = m.cant_columnas();
       
        _m_filas = 2 * m.cant_filas() *  m.cant_columnas();
        _m_colum = (m.cant_filas() + 1) * (m.cant_columnas() + 1) - 1;
       
        for(unsigned int i = 0; i < _m_filas; i++)
                _m.push_back(vpid());

        unsigned int fila_i_x = 0;
        unsigned int colu_i_x = 0;
       
        unsigned int fila_i_y = 1;
        unsigned int colu_i_y = 0;

        for(unsigned int i = 0; i < m.cant_filas(); i++)
        {
                for(unsigned int j = 0; j < m.cant_columnas(); j++)
                {
                        //aca los elementos correspondientes a las filas x
                        _m_set_elem(fila_i_x, colu_i_x, -m[i][j].trd());
                        _m_set_elem(fila_i_x, colu_i_x +  m.cant_columnas() + 1, m[i][j].trd());
                                               
                        fila_i_x += 2;
                        colu_i_x++;
       
                        //aca los elementos correspondientes a las filas y
                        _m_set_elem(fila_i_y, colu_i_y, -m[i][j].trd());
                        _m_set_elem(fila_i_y, colu_i_y + 1, m[i][j].trd());
                               
                        fila_i_y += 2;
                        colu_i_y++;
                       
                        //si termine de recorrer una fila en la matriz original debo avanzar una columna mas   
                        if(j ==  m.cant_columnas() - 1)
                        {
                                colu_i_x++;
                                colu_i_y++;
                        }                      
                }
        }
}

/* IMPORTANTE:
 * Esta construccion podria acelerarse si aprovechamos el hecho de que la matriz resultante es simetrica */

void matriz_esp::_construir_mtxm()
{
        /* La matriz resultado tiene la cantidad de filas de Mt por cantidad de columnas de M, que en ambos casos es la cantidad de columnas de M */
       
        _mtxm_filas = _mtxm_columnas = _m_colum;
       
        for(unsigned int i = 0; i < _m_colum; i++)
        {
                _mtxm_pf.push_back(vpid());
                _mtxm_pc.push_back(vpid());
        }
               
        /* Sabemos que la matriz m tiene 2 elementos por cada fila (sino esto no funciona!) */
       
        for(unsigned int i = 0; i < _m_filas; i++)
        {
                /* Estoy en la iesima fila, tengo que hacer 4 multiplicaciones */
               
                vpid fila = _m[i];
                assert(fila.size() == 2);
               
                pair<unsigned int, double> prim = fila[0];
                pair<unsigned int, double> segu = fila[1];
                unsigned int pos_prim = prim.first;
                unsigned int pos_segu = segu.first;
               
                _mtxm_add_elem(pos_prim, pos_prim, prim.second * prim.second);
                _mtxm_add_elem(pos_prim, pos_segu, prim.second * segu.second);
                _mtxm_add_elem(pos_segu, pos_prim, segu.second * prim.second);
                _mtxm_add_elem(pos_segu, pos_segu, segu.second * segu.second);           
        }
}

void matriz_esp::_construir_b(matriz<tripla_dbl> & m)
{
        vector<double> b_temp;
       
        for(unsigned int i = 0; i < m.cant_filas(); i++)
        {
                for(unsigned int j = 0; j < m.cant_columnas(); j++)
                {
                        b_temp.push_back(-(m[i][j]).fst());
                        b_temp.push_back(-(m[i][j]).snd());
                }      
        }
               
        //cout << "b_temp.size()   " << b_temp.size() << endl;
       
        b = mult_mt(b_temp);                            // Multiplica Mt * B
       
        //cout << "b.size()        " << b.size() << endl;
       
        x = vector<double>(b.size(), 0);        // Inicializa el vector de incognitas  
}

/** </METODOS PRIVADOS> **/

/** </CONSTRUCTORES> **/



/** <METODOS PUBLICOS DE M> **/

/* Multiplica Mt por el vector V */

vector<double> matriz_esp::mult_mt(vector<double> &v)
{
        if(_m_filas != v.size())
                cout << "mult_mt(): Las dimensiones no concuerdan!" << endl;
       
        vector<double> res(_m_colum, 0);
       
        for(unsigned int i = 0; i < _m_filas; i++)
        {
                vpid fila = _m[i];
                assert(fila.size() == 2);
               
                pair<unsigned int, double> prim = fila[0];
                pair<unsigned int, double> segu = fila[1];
                unsigned int pos_prim = prim.first;
                unsigned int pos_segu = segu.first;
               
                res[pos_prim] += v[i] * prim.second;
                res[pos_segu] += v[i] * segu.second;   
        }
       
        return res;    
}

/* Obtiene un elemento de la matriz M */

double matriz_esp::get_elem_m(unsigned int fila, unsigned int columna)
{
        return _m_get_elem(fila, columna);
}

/* Muestra la matriz M por pantalla */

void matriz_esp::mostrar_m()
{
        for(unsigned int i = 0; i < _m_filas; i++)
        {
                for(unsigned int j = 0; j < _m_colum; j++)
                        cout << _m_get_elem(i, j) << " ";
                       
                cout << endl;  
        }
}

void matriz_esp::mostrar_b()
{
        for(unsigned int i = 0; i < b.size(); i++)
                cout << b[i] << " ";
               
        cout << endl;
}

/** </METODOS PUBLICOS DE M> **/



/** <METODOS PRIVADOS DE M> **/

double matriz_esp::_m_get_elem(unsigned int fila, unsigned int columna)
{
        /* NOTA:
         * Este codigo podria cambiarse porque _m[fila].size() deberia ser siempre 2,
         * pero como funciona igual de rapido, lo dejamos asi */
       
        for(unsigned int i = 0; i < _m[fila].size(); i++)
                if(_m[fila][i].first == columna)
                        return _m[fila][i].second;
       
        /* Si el elemento no esta, significa que es igual a 0 */
       
        return 0;
}

/* Seteo un elemento en M, que sabemos va a tener 2 elementos por fila */

bool matriz_esp::_m_set_elem(unsigned int fila, unsigned int columna, double valor)
{
        if((_m_filas - 1 < fila) || (_m_colum - 1 < columna))
        {
                cout << "_m_set_elem: La posicion (" << fila << ", " << columna << ") esta fuera de rango en la matriz M!" << endl;
                return false;
        }
       
        if(esCero(valor))
        {
                cout << "_m_set_elem: La matriz M no deberia setear ceros!" << endl;
                return false;
        }
       
        if(_m[fila].size() == 2)        /* ERROR */
        {
                if((_m[fila][0].first == columna) || (_m[fila][1].first == columna))
                        cout << "_m_set_elem: La posicion (" << fila << ", " << columna << ") ya existe en la matriz M y no deberia sobreescribirse! " << endl;
                else
                        cout << "_m_set_elem: La matriz ya tiene dos elementos en la fila " << fila << "!" << endl;
                       
                return false;
        }
                       
        /* Si sigo por aca es porque esta todo bien */
       
        if(_m[fila].size() == 1)
        {
                if(_m[fila][0].first < columna) _m[fila].push_back(pair<unsigned int, double>(columna, valor));
                else
                {
                        if(_m[fila][0].first == columna)
                        {
                                cout << "_m_set_elem: La posicion (" << fila << ", " << columna << ") ya existe en la matriz M y no deberia sobreescribirse! " << endl;
                                return false;
                        }
                       
                        /* Si sigo por aca es xq tengo que reordenar */
                       
                        pair<unsigned int, double> temp = _m[fila][0];
                        _m[fila][0] = pair<unsigned int, double>(columna, valor);
                        _m[fila][1] = temp;
                }
        }
       
        if(_m[fila].size() == 0)
                _m[fila].push_back(pair<unsigned int, double>(columna, valor));

        return true;   
}

/** </METODOS PRIVADOS DE M> **/



/** <METODOS PUBLICOS DE MtxM> **/

double matriz_esp::get_elem_mtxm(unsigned int fila, unsigned int columna)
{
         return _mtxm_get_elem_pf(fila, columna);
}

/* Muestra la matriz MtxM como le gusta a Martin */

void matriz_esp::mostrar_mtxm()
{
        for(unsigned int i = 0; i < _mtxm_filas ; i++)
        {
                for(unsigned int j = 0; j < _mtxm_columnas; j++)
                        cout << _mtxm_get_elem_pf(i, j) << " ";
                           
                cout << endl;
        }
       
}

void matriz_esp::mostrar_mtxm_pf()
{
        for(unsigned int i = 0; i < _m_colum ; i++)
        {
                for(unsigned int j = 0; j < _m_colum ; j++)
                {
                        cout << setw(12);
                        cout << setprecision(10);
                        cout << _mtxm_get_elem_pf(i,j) << " " ;
                }
                       
                cout << endl;
        }
}

void matriz_esp::mostrar_mtxm_pc()
{
        for(unsigned int i = 0; i < _m_colum ; i++)
        {
                for(unsigned int j = 0; j < _m_colum ; j++)
                {
                        cout << setw(12);
                        cout << setprecision(10);
                        cout << _mtxm_get_elem_pc(i,j) << " " ;
                }
                       
                cout << endl;
        }
}

/* Resuelve el sistema de ecuaciones y da el resultado */

vector<double> matriz_esp::resolver()
{
        _mtxm_triangular();
                       
        //cout << endl;
        //mostrar_mtxm_pc();
        //cout << endl;
        //mostrar_mtxm_pf();
        //cout << endl;
        //mostrar_b();
       
        _mtxm_backsust();
                                       
        return x;
}


/** </METODOS PUBLICOS DE MtxM> **/



/** <METODOS PRIVADOS DE MtxM> **/

/** ---------------------------- GETS ------------------------------ **/

double matriz_esp::_mtxm_get_elem_pf(unsigned int fila, unsigned int columna)
{
        /* NOTA:
         * Aca no puedo hacer mejor esto porque a priori no sabemos la
         * cantidad de elementos por fila ni columna */
         
        for(unsigned int i = 0; i < _mtxm_pf[fila].size(); i++)
                if(_mtxm_pf[fila][i].first == columna)
                        return _mtxm_pf[fila][i].second;
       
        /* Si el elemento no esta, significa que es igual a 0 */
       
        return 0;
}

double matriz_esp::_mtxm_get_elem_pc(unsigned int fila, unsigned int columna)
{
        for(unsigned int i = 0; i < _mtxm_pc[columna].size() ; i++)
                if(_mtxm_pc[columna][i].first == fila)
                        return _mtxm_pc[columna][i].second;
                               
        return 0;
}

/** ---------------------------- SETS ------------------------------ **/

bool matriz_esp::_mtxm_set_elem(unsigned int fila, unsigned int columna, double valor)
{
        if((_m_colum - 1 < fila) || (_m_colum - 1 < columna))           /* La cantidad de columnas de M es la cantidad de filas y columnas de MtxM */
        {
                cout << "_mtxm_set_elem: La posicion (" << fila << ", " << columna << ") esta fuera de rango en la matriz M!" << endl;
                return false;
        }
       
        /* IMPORTANTE:
         * Hago todo aca sin llamar a funciones auxiliares para evitar recorrer la misma estructura varias veces */
       
        /* Primero me paro en la posicion que quiero setear */
       
        /* Agregamos por filas */
       
        _mtxm_set_elem_pf(fila, columna, valor);
        _mtxm_set_elem_pc(fila, columna, valor);
               
        return true;
}

bool matriz_esp::_mtxm_set_elem_pf(unsigned int fila, unsigned int columna, double valor)
{
        for(unsigned int i = 0; i < _mtxm_pf[fila].size(); i++)
        {
                if(_mtxm_pf[fila][i].first == columna)
                {
                        /* Si encontre la posicion que quiero, veo que tengo que hacer... */
                       
                        if(esCero(valor))
                                _mtxm_pf[fila].erase(_mtxm_pf[fila].begin() + i);
                        else
                                _mtxm_pf[fila][i].second = valor;
                       
                        return true;                   
                }
                if(_mtxm_pf[fila][i].first > columna) break;    /* Esto funciona porque las columnas estan ordenadas */
        }
               
        /* Si sigo por aca es xq la posicion no existe */
       
        if(esCero(valor)) return true;  /* Si lo que queria hacer era borrar, no hay nada que hacer */
       
        _mtxm_pf[fila].push_back(pair<unsigned int, double>(columna, valor));
        sort(_mtxm_pf[fila].begin(), _mtxm_pf[fila].end(), primeraComponenteMenor);     /* Ordeno el vector */
        return true;           
}

bool matriz_esp::_mtxm_set_elem_pc(unsigned int fila, unsigned int columna, double valor)
{
        for(unsigned int i = 0; i < _mtxm_pc[columna].size(); i++)
        {
                if(_mtxm_pc[columna][i].first == fila)
                {
                        /* Si encontre la posicion que quiero, veo que tengo que hacer... */
                       
                        if(esCero(valor))
                                _mtxm_pc[columna].erase(_mtxm_pc[columna].begin() + i);
                        else
                                _mtxm_pc[columna][i].second = valor;
                       
                        return true;                   
                }
                if(_mtxm_pc[columna][i].first > fila) break;    /* Esto funciona porque las filas estan ordenadas */
        }
               
        /* Si sigo por aca es xq la posicion no existe */
       
        if(esCero(valor)) return true;  /* Si lo que queria hacer era borrar, no hay nada que hacer */
       
        _mtxm_pc[columna].push_back(pair<unsigned int, double>(fila, valor));
        sort(_mtxm_pc[columna].begin(), _mtxm_pc[columna].end(), primeraComponenteMenor);       /* Ordeno el vector */
        return true;           
}

/** ---------------------------- ADDS ------------------------------ **/

/* El 'add' de esta funcion es de "sumar" no de "añadir" */

bool matriz_esp::_mtxm_add_elem(unsigned int fila, unsigned int columna, double valor)
{
        if((_m_colum - 1 < fila) || (_m_colum - 1 < columna))
        {
                cout << "_mtxm_add_elem: La posicion (" << fila << ", " << columna << ") no existe en la matriz MtxM!" << endl;
                return false;
        }
       
        if(esCero(valor)) return true; /* Si esta sumando cero, no hago nada */
       
        /* IMPORTANTE:
         * Hago todo aca sin llamara funciones auxiliares para evitar recorrer la misma estructura varias veces.
         * Si, queda el codigo duplicado (o hasta triplicado) pero ahorro tiempo. Y el tiempo es oro... */
       
        _mtxm_add_elem_pf(fila, columna, valor);
        _mtxm_add_elem_pc(fila, columna, valor);
       
        return true;
}

bool matriz_esp::_mtxm_add_elem_pf(unsigned int fila, unsigned int columna, double valor)
{
        for(unsigned int i = 0; i < _mtxm_pf[fila].size(); i++)
        {
                if(_mtxm_pf[fila][i].first == columna)
                {
                        double nuevo_valor = _mtxm_pf[fila][i].second + valor;
                        if(esCero(nuevo_valor))
                               
                                _mtxm_pf[fila].erase(_mtxm_pf[fila].begin() + i);
                       
                        else
                                _mtxm_pf[fila][i].second = nuevo_valor;

                        return true;                   
                }
                if(_mtxm_pf[fila][i].first > columna) break;    /* Esto funciona porque las columnas estan ordenadas */
        }
       
        /* Si sigo por aca es porque la posicion no existe, por lo que vale 0 */
       
        _mtxm_pf[fila].push_back(pair<unsigned int, double>(columna, valor));
        sort(_mtxm_pf[fila].begin(), _mtxm_pf[fila].end(), primeraComponenteMenor);                     /* Ordeno el vector */
               
        return true;
}

bool matriz_esp::_mtxm_add_elem_pc(unsigned int fila, unsigned int columna, double valor)
{
        for(unsigned int i = 0; i < _mtxm_pc[columna].size(); i++)
        {
                if(_mtxm_pc[columna][i].first == fila)
                {
                        double nuevo_valor = _mtxm_pc[columna][i].second + valor;
                        if(esCero(nuevo_valor))
                                _mtxm_pc[columna].erase(_mtxm_pc[columna].begin() + i);                
                        else
                                _mtxm_pc[columna][i].second = nuevo_valor;
                               
                        return true;                   
                }
                if(_mtxm_pc[columna][i].first > fila) break;    /* Esto funciona porque las columnas estan ordenadas */
        }
       
        /* Si sigo por aca es porque la posicion no existe, por lo que vale 0 */
       
        _mtxm_pc[columna].push_back(pair<unsigned int, double>(fila, valor));
        sort(_mtxm_pc[columna].begin(), _mtxm_pc[columna].end(), primeraComponenteMenor);                       /* Ordeno el vector */
       
        return true;
}

/** ---------------------------- DELS ------------------------------ **/

bool matriz_esp::_mtxm_del_elem(unsigned int fila, unsigned int columna)
{
        if((_m_colum - 1 < fila) || (_m_colum - 1 < columna))
        {
                cout << "_mtxm_del_elem: La posicion (" << fila << ", " << columna << ") no existe en la matriz MtxM!" << endl;
                return false;
        }
       
        return _mtxm_set_elem(fila, columna, 0);
}

/** ---------------------- OTRAS FUNCIONES ------------------------- **/

/* Funcion que triangula la matriz MtxM */

void matriz_esp::_mtxm_triangular()
{
        if(mostrar_debug)
        {
                cout << "Triangulando matriz esparsa..." << endl << endl;
                cout << "  0 %" << endl;
        }
       
        int porcentaje_ant = 0;
       
        /* Version mas copada de gauss by Martin C. y  Martin S. */
               
        for(unsigned int columna_a_anular = 0; columna_a_anular < _m_colum; columna_a_anular++)
        {               //cout << "columna " << columna_a_anular << endl;
                if(mostrar_debug)       // Mensaje de debug, no darle bola
                {
                        double porcentaje = (columna_a_anular * 100) / _m_colum;
                        if(porcentaje > porcentaje_ant)
                        {
                                porcentaje_ant = porcentaje;
                                cout << setw(3);
                                cout << porcentaje << " %" << endl;
                        }
                }              
               
                /* Pongo en 0 la columna a anular y obtengo los multiplicadores */
               
                double elem_diagonal = 0;
               
                for(unsigned int i = 0; i < _mtxm_pc[columna_a_anular].size(); i++)
                {
                        if(_mtxm_pc[columna_a_anular][i].first == columna_a_anular)
                        {
                                elem_diagonal = _mtxm_pc[columna_a_anular][i].second;
                                break;
                        }
                       
                        if(_mtxm_pc[columna_a_anular][i].first > columna_a_anular)
                                break;
                }              
                       
                if(esCero(elem_diagonal))
                {
                        cout << "HAY UN CERO EN LA DIAGONAL " << columna_a_anular << " " << columna_a_anular << endl;
                }
                else
                {      
                        /* El elemento de la diagonal es distinto de 0, creamos los multiplicadores */
                       
                        vector<pid> multiplicadores;

                        for(unsigned int i = _mtxm_pc[columna_a_anular].size() - 1; i > 0; i--)
                        {
                                unsigned int fila = _mtxm_pc[columna_a_anular][i].first;
                               
                                if(fila <= columna_a_anular) break;             /* Si estoy en la diagonal o por arriba, salgo */
                               
                                /* Si sigo por aca es xq estoy por abajo de la diagonal */
                                                               
                                multiplicadores.push_back(pair<unsigned int, double>(fila, (-(_mtxm_pc[columna_a_anular][i].second) / elem_diagonal)));
                               
                                _mtxm_pc[columna_a_anular].pop_back(); // borro el ultimo porque queremos poner ceros en la columna                    
                        }

                        /* Ya tenemos los multipicadores, hay que usarlos. Iteramos las columnas aplicando los multiplicadores */
                       
                        for(unsigned int columna_actual = columna_a_anular + 1; columna_actual < _m_colum ; columna_actual++)
                        {
                                /* Queremos el elemento que esta en la misma fila que el elemento de la diagonal, pero en la columna actual */
                               
                                double constante_multiplicador = 0;
                                                       
                                for(unsigned int i = 0; i < _mtxm_pc[columna_actual].size(); i++)
                                {
                                        if(_mtxm_pc[columna_actual][i].first == columna_a_anular)       //  Encontre el elemento que esta en la columna a modificar y en la fila de la diagonal
                                        {
                                                constante_multiplicador = _mtxm_pc[columna_actual][i].second;
                                                break;
                                        }
                                       
                                        if(_mtxm_pc[columna_actual][i].first > columna_a_anular)        // Me pase y no lo encontre
                                                break;
                                }
                               
                                // Si la constante del multiplicador es 0, significa que las posiciones de la columna actual que estan por debajo de la diagonal no se modifican
                               
                                if(constante_multiplicador == 0)
                                        continue;
                       
                                // Si llegamos hasta aca tenemos que modificar los valores que estan por debajo del elemento de la diagonal en la columna actual
                                                               
                                unsigned int puntero_a_columna = _mtxm_pc[columna_actual].size() - 1;                          
                               
                                for(unsigned int puntero_multiplicadores = 0; puntero_multiplicadores < multiplicadores.size(); puntero_multiplicadores++)
                                {
                                        // El vector de multiplicadores se recorre en orden y la columna al reves
                                       
                                        unsigned int fila_multi = multiplicadores[puntero_multiplicadores].first;
                                        double elemento_actual = 0;                                    
                                       
                                        while(true)
                                        {
                                                if(_mtxm_pc[columna_actual][puntero_a_columna].first == fila_multi)
                                                {
                                                        elemento_actual = _mtxm_pc[columna_actual][puntero_a_columna].second;
                                                        break;                                                 
                                                }
                                               
                                                if(_mtxm_pc[columna_actual][puntero_a_columna].first < fila_multi)
                                                        break;

                                                puntero_a_columna--;
                                        }
                                       
                                        double valor_a_poner = elemento_actual + constante_multiplicador * multiplicadores[puntero_multiplicadores].second;
                                                                               
                                        if(elemento_actual == 0)
                                        {
                                                if(puntero_a_columna == (_mtxm_pc[columna_actual].size() - 1))
                                                {
                                                        _mtxm_pc[columna_actual].push_back(pair<unsigned int,double>(fila_multi,valor_a_poner));
                                                }
                                                else
                                                {
                                                        vector< pid>::iterator it;
                                                        it = _mtxm_pc[columna_actual].begin();
                                                        it = _mtxm_pc[columna_actual].insert(it + puntero_a_columna + 1, pair<unsigned int,double>(fila_multi, valor_a_poner));
                                                }
                                                continue;
                                        }
                                       
                                        if(esCero(valor_a_poner))
                                        {
                                                _mtxm_pc[columna_actual].erase(_mtxm_pc[columna_actual].begin()+ puntero_a_columna );
                                                continue;
                                        }
                                       
                                        /* Si esta todo bien, solo modificamos el valor que existia */
                                                               
                                        _mtxm_pc[columna_actual][puntero_a_columna].second = valor_a_poner;                                    
                                }                                      
                        }      
                       
                        /* Ya modificamos todas las columnas, ahora modificamos el termino independiente */
                       
                        double constante_multiplicador_b = b[columna_a_anular];         // Recordar que columna a anular es la fila donde esta el elemento de la diagonal
                                                                               
                        // Al igual que antes, si la constante del multiplicador es 0, no se cambia nada
                       
                        if(esCero(constante_multiplicador_b))
                                continue;
               
                        // Recorro el vector de multiplicadores, y modifico las posiciones relacionadas con el multiplicador
                       
                        for(unsigned int puntero_multiplicadores = 0; puntero_multiplicadores < multiplicadores.size(); puntero_multiplicadores++)
                        {
                                unsigned int fila_multi = multiplicadores[puntero_multiplicadores].first;                              
                                b[fila_multi] += constante_multiplicador_b * multiplicadores[puntero_multiplicadores].second;
                        }
                }
        }
       
        if(mostrar_debug)
        {
                cout << "100 %" << endl << endl;
                cout << "Matriz esparsa triangulada." << endl;
        }
       
        _mtxm_reestablecer_inv_pf();
}

/* Despeja las incognitas con backward sustitution */

void matriz_esp::_mtxm_backsust()
{
        if(mostrar_debug)
        {
                cout << "Despejando incognitas..." << endl << endl;
                cout << "  0 %" << endl;
        }
       
        int porcentaje_ant = 0;
       
        for(int i = _mtxm_filas - 1; i >= 0; i--)
        {
                if(mostrar_debug)       // Mensaje de debug, no darle bola
                {
                        double porcentaje = ((_mtxm_filas - i - 1) * 100) / _mtxm_filas;
                        if(porcentaje > porcentaje_ant)
                        {
                                porcentaje_ant = porcentaje;
                                cout << setw(3);
                                cout << porcentaje << " %" << endl;
                        }
                }
               
                x[i] = b[i];
               
                int it;
                unsigned int col;              
               
                for(it = _mtxm_pf[i].size() - 1; it >= 0; it--)         // Recorro la fila de atras para adelante
                {                      
                        // Agarro el elemento en la posicion it
                        // Veo en que columna esta el elemento
                        // - Si la columna es mayor que i, procedo
                        // - Si la columna es menor o igual que i, corto
                       
                        if((col = _mtxm_pf[i][it].first) > i)
                                x[i] -= _mtxm_pf[i][it].second * x[col];
                               
                        else
                        {
                                // cout << i << " " << _mtxm_pf[i][it].first << endl;   // Para ver que efectivamente paramos en la diagonal
                                break;
                        }
                }
               
                // Como MtxM no tiene ceros en la diagonal, en _mtxm_pf[i][it] tengo el elemento de la diagonal
                               
                x[i] /= _mtxm_pf[i][it].second;
        }
       
        /* Despeje clasico

    for(int i = _mtxm_filas - 1; i >= 0; i--)
    {
                x[i] = b[i];
                for(int j = i + 1; j < (int) _mtxm_columnas; j++)
                        x[i] -= _mtxm_get_elem_pc(i, j) * x[j];
               
                x[i] /= _mtxm_get_elem_pc(i, i);
        }
        */
       
        if(mostrar_debug)
        {
                cout << "100 %" << endl << endl;
                cout << "Incognitas despejadas." << endl;
        }
}

/* Esta funcion hace que _mtxm_pf sea consistente con _mtxm_pc */

void matriz_esp::_mtxm_reestablecer_inv_pf()
{
        if(mostrar_debug)
        {
                cout << "Reestableciendo invariante..." << endl << endl;
                cout << "  0 %" << endl;
        }
       
        int porcentaje_ant = 0;
       
        // Primero limpio la estructura
       
        for(unsigned int i = 0; i < _mtxm_pf.size(); _mtxm_pf[i].clear(), i++)
                ;
               
        // Ahora la completo
       
        for(unsigned int col = 0; col < _mtxm_pc.size(); col++)
        {
                if(mostrar_debug)       // Mensaje de debug, no darle bola
                {
                        double porcentaje = (col * 100) / _mtxm_pc.size();
                        if(porcentaje > porcentaje_ant)
                        {
                                porcentaje_ant = porcentaje;
                                cout << setw(3);
                                cout << porcentaje << " %" << endl;
                        }
                }
               
                for(unsigned int i = 0; i < _mtxm_pc[col].size(); i++)
                {
                        unsigned int fila = _mtxm_pc[col][i].first;
                        double elem = _mtxm_pc[col][i].second;                 
                        _mtxm_pf[fila].push_back(pid(col, elem));
                }
        }
       
        if(mostrar_debug)
        {
                cout << "100 %" << endl << endl;
                cout << "Invariante reestablecido." << endl;
        }
}


/** </METODOS PRIVADOS DE MtxM> **/
